function [rationum ratioMean ratioStd]=image_contrast(im_pre,im_after,xpeak,ypeak,radial,bgmask)
%if rationum and ratioMean more than 1,pictures become better;if ratioStd
%less than 1,pictures become better.


min_value=30;
max_value=900;
diff_extend=10;
extend_inside=2;
xpeak_length=length(xpeak);
rationum=zeros(xpeak_length,1);
ratioMean=zeros(xpeak_length,1);

%change number of row and column
im_size=size(im_pre);
imprec=reshape(im_pre,im_size(1)*im_size(2),1);
imafterc=reshape(im_after,im_size(1)*im_size(2),1);
k = find(bgmask(:)==1);  % Index k specifying background region

im_prebgmask=imprec(k,1);
im_afterbgmask=imafterc(k,1);
% im_prebgmask=im_prebgmask(:);
% im_afterbgmask=im_afterbgmask(:);
MatrixPreMeanUnROI=mean(im_prebgmask);
MatrixAfterMeanUnROI=mean(im_afterbgmask);
MatrixPreStdUnROI=std(im_prebgmask);
MatrixAfterStdUnROI=std(im_afterbgmask);
ratioStd=MatrixAfterStdUnROI/MatrixPreStdUnROI;

for i=1:xpeak_length
    if radial(i,1)==0
       break; 
    end
    xLeftLimit=xpeak(i,1)-radial(i,1);
    xRightLimit=xpeak(i,1)+radial(i,1);
    yBottomLimit=ypeak(i,1)-radial(i,1);
    yTopLimit=ypeak(i,1)+radial(i,1);
    
    %limit range of x,y
    if xLeftLimit-diff_extend<min_value;xLeftLimit=diff_extend+min_value;end
    if xRightLimit+diff_extend>max_value;xRightLimit=max_value-diff_extend;end
    if yTopLimit-diff_extend<min_value;yTopLimit=min_value+diff_extend;end
    if yBottomLimit+diff_extend>max_value;yBottomLimit=max_value-diff_extend;end
    
    Matrix_xLeftPre=im_pre(yBottomLimit:yTopLimit,xLeftLimit-diff_extend:xLeftLimit+radial/extend_inside);
    Matrix_xRightPre=im_pre(yBottomLimit:yTopLimit,xRightLimit-radial/extend_inside:xRightLimit+diff_extend);
    Matrix_yTopPre=im_pre(yBottomLimit-diff_extend:yBottomLimit+radial/extend_inside,xLeftLimit:xRightLimit);
    Matrix_yBottomPre=im_pre(yTopLimit-radial/extend_inside:yTopLimit+radial/diff_extend,xLeftLimit:xRightLimit);
    
    Matrix_xLeftAfter=im_after(yBottomLimit:yTopLimit,xLeftLimit-diff_extend:xLeftLimit+radial/extend_inside);
    Matrix_xRightAfter=im_after(yBottomLimit:yTopLimit,xRightLimit-radial/extend_inside:xRightLimit+diff_extend);
    Matrix_yTopAfter=im_after(yBottomLimit-diff_extend:yBottomLimit+radial/extend_inside,xLeftLimit:xRightLimit);
    Matrix_yBottomAfter=im_after(yTopLimit-radial/extend_inside:yTopLimit+radial/diff_extend,xLeftLimit:xRightLimit);
    
    MatrixPreMeanROI=mean(mean(im_pre(yBottomLimit:yTopLimit,xLeftLimit:xRightLimit)));
    MatrixAfterMeanROI=mean(mean(im_after(yBottomLimit:yTopLimit,xLeftLimit:xRightLimit)));
    
    ratioMeanPre=MatrixPreMeanROI-MatrixPreMeanUnROI;
    ratioMeanAfter=MatrixAfterMeanROI-MatrixAfterMeanUnROI;
    
    ratioMean(i,1)=ratioMeanAfter/ratioMeanPre;
    
    
    diff_xLeftPre=diff((Matrix_xLeftPre).');
    diff_xRightPre=diff((Matrix_xRightPre).');
    diff_yTopPre=diff(Matrix_yTopPre);
    diff_yBottomPre=diff(Matrix_yBottomPre);
    
    diff_xLeftAfter=diff((Matrix_xLeftAfter).');
    diff_xRightAfter=diff((Matrix_xRightAfter).');
    diff_yTopAfter=diff(Matrix_yTopAfter);
    diff_yBottomAfter=diff(Matrix_yBottomAfter);
    
    diff_xLeftPreMean=mean(mean(diff_xLeftPre));
    diff_xRightPreMean=mean(mean(diff_xRightPre));
    diff_yTopPreMean=mean(mean(diff_yTopPre));
    diff_yBottomPreMean=mean(mean(diff_yBottomPre));
    
    %number of >differential mean
    num_xLeftPre=length(find(diff_xLeftPre>diff_xLeftPreMean));
    num_xRightPre=length(find(diff_xRightPre>diff_xRightPreMean));
    num_yTopPre=length(find(diff_yTopPre>diff_yTopPreMean));
    num_yBottomPre=length(find(diff_yBottomPre>diff_yBottomPreMean));
    
    num_xLeftAfter=length(find(diff_xLeftAfter>diff_xLeftPreMean));
    num_xRightAfter=length(find(diff_xRightAfter>diff_xRightPreMean));
    num_yTopAfter=length(find(diff_yTopAfter>diff_yTopPreMean));
    num_yBottomAfter=length(find(diff_yBottomAfter>diff_yBottomPreMean));
    
    numPre=num_xLeftPre+num_xRightPre+num_yTopPre+num_yBottomPre;
    numAfter=num_xLeftAfter+num_xRightAfter+num_yTopAfter+num_yBottomAfter;
    rationum(i,1)=numAfter/numPre;
end






